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Abstract 

We model two-dimensional systems of granular aggregates con- 
fined between two planes and demonstrate that at a critical grain 
volume fraction an abrupt rigidity transition occurs. This transition 
is observed both in static and shear tests. The grain volume fraction 
at which the transition occurs, u c , decreases with increasing friction 
between the grains. Densely packed grains, with a volume fraction 
v > u c , display an elastic-plastic rheology. Dilute packings, with 
v < u c , display gas- like characteristics. It is shown that when vol- 
ume fraction is allowed to change freely (using constant normal stress 
boundary condition), it evolves spontaneously to v c under a wide range 
of boundary conditions, exhibiting 'self-organized criticality'. 

Granular media is a fundamental, yet not well understood, complex sys- 
tem with wide ranging applications to technological and natural systems. In 
recent years there has been much work on granular dynamics, with emphasis 
on how behavior of grain aggregates may resemble solids, liquids or gases 
[|TJ. Already Reynolds in 1885 H noted that loosely packed sands deform 
easily as fluids while dense packings resist shear as solids. The two phases 
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are traditionally treated separately by 'kinetic-gas' approaches for loosely 
packed grains [|J and elasto-plastic (often using associated plasticity f|) the- 
ories for dense soils. In this Letter we numerically investigate the transition 
between solid and gas behaviors in granular aggregates, and show that the 
phase boundary between gas and solid is also an attracting point, to which 
systems naturally evolve. We numerically model grain aggregates using a 
version of the popular 'discrete element method' || which treats grains as 
inelastic disks with rotational and translational degrees of freedom. Two 
grains of radius Ri and Rj undergo an inelastic interaction when the dis- 
tance separating them m is less than the sum of their radii. During the 
interaction the ith grain experiences a contact force that has both shear and 
normal components: 

Fij(t) = [k n (Ri + Rj - Tij) - 7m(rij ■ n)]n + [min(k s As, fi(F ■ n))]s (1) 

where n = (r^ ■ x, ■ y)/ry, § = (r^ ■ y, — ■ x)/ry, are the unit vectors 
in the normal and tangential directions respectively. k n , k s are the normal 
and shear elastic constants, m is the grain mass, 7 is a damping coefficient 
ensuring inelasticity of the interaction, fj, is the surface friction coefficient, 
and As is the shear displacement since the initial contact of two particles. 
Force is integrated through time to calculate grain position and velocity. 
For collisions governed by equation (|l|) energy loss is governed by a normal 
restitution coefficient e n = exp(— 7^/2), where t co i = ir(2k n /m — 7 2 /4)~ 1 / 2 
is the collision time, and by frictional work which depends on the amount 
of real slip (after the elastic limit is reached) and the frictional shear force, 
fj,F n . Time is measured in units of undissipated elastic wave travel time 
to = \Jm/k n , and distance in units of average disk diameter xq = 2R. In 
simulations presented here k n = l,j = l,k a = 0.5 and m = 1. 

In this Letter we investigate three related problems. The first problem 
is the general characteristics of static granular aggregates after compaction 
between two parallel plates. The second and third problems are the behavior 
of the same confined aggregates during shear using two different boundary 
conditions: constant distance between shearing walls, termed constant vol- 
ume boundary conditions (CVBC), and constant normal stress applied to 
the wall, termed constant force boundary conditions (CFBC). Simulations 
are performed on square systems with n disks. The top and bottom edges of 
the box are composed of grains glued together to form rigid rough walls of 
length I (Figure [I]). The box is periodic in the horizontal direction. Grain 
radii are randomly drawn from a Gaussian distribution that peaks at R, with 
a standard deviation of 0.5/2. Polydispersivity is introduced to discourage 
ordering effects. The system is initiated as tall loosely packed box, which is 
compacted vertically by normal stresses to a height I. After compaction the 
horizontal walls are allowed to move freely in the horizontal direction to relax 
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forces. However, unrelaxed normal forces can be maintained on the walls in 
this case, since we do not allow the walls to move vertically, using CVBC. 
Global rearrangements during compaction and relaxation ensure (local) min- 
imal energy configuration, as would occur during natural compaction. 

After compaction and relaxation we measure properties of static configu- 
rations at different solid volume fractions, v = J2i=i ^R^/l 2 , ranging between 
0.75 to 0.96. Measurements of three parameters, (a) the average number of 
grains touching (i.e. exerting a force on) a grain in the interior of the box, 
termed the coordination number Z, (b) The normal stress (normal force per 
unit length) N operating on the upper and lower confining walls, and (c) the 
systems shear modulus, G, are presented in Figure |2|. All measures show an 
abrupt change in behavior at a critical volume fraction z/ c , which depends 
on the value of friction prescribed between the grains, but not on the sys- 
tem size. The coordination number is approximately zero for v < v c . At 
v c the coordination number abruptly jumps and then increases as an em- 
pirical power-law, Z oc [y — v cp ) a ^ , where the subscript fi denotes different 
friction values used in simulations. Power-law fits (in solid curves) yields 
z/ c0 = 0.83 ± 0.01, a = 0.5 ± 0.1 and u c0 5 = 0.805 ± 0.01, a . 5 = 0.3 ± 0.05. 
Critical behavior with values of z/ c0 = 0.82 ± 0.02 were obtained for hard 
smooth disks ||, and visco-elastic 2D bubbles |7[ for mono-dispersed and 
polydispersed systems, demonstrating that v c is fairly independent of the 
disk size distribution, and the interaction law in absence of friction. The 
difference in v c and a between frictional and smooth grains occurs because 
frictional grains tend to 'stick', and thus cannot achieve the lower energy 
configuration of smooth disks. 

The fact that N and Z (Figure |2|a,b) are zero for v < u c , indicates that 
below v c grains do not touch. At v c grains first touch and elastically repel each 
other, exerting normal stresses on the walls. For v > u c , normal stress follows 
iV oc Z[y — v c ) (shown in solid lines), a slightly modified form (accounting for 
the phase transition) of predictions from standard models of densely packed 
elastic disks || . The observed transition is identified as a macroscopic rigidity 
transition in Figure |2]c. The system's elastic shear modulus G is obtained 
by imposing a small homogeneous shear step strain of a magnitude e and 
measuring the resulting shear stress on the walls o [G = o jt is independent 
of e for e < 10~ 4 of the grain radius, here we use e = 10 -5 ). The procedure 
follows that outlined in |7]] . Figure [|c shows G increasing from zero to a finite 
value, as the system passes through the phase transition at v = v c . Physically, 
the rigidity should depend on average number of springs per disk, yet not 
on their compression, i.e. G oc Z (solid curves). The rigidity-transition thus 
identified is a result of global geometrical constraints of grain packings, with 
the coordination number identified as the order parameter for the transition. 
It was previously shown that v = 0.83 ± 0.02 marks the upper limit of 
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compacity of disordered packings of smooth hard mono-sized and poly-sized 
disks ||, where long-range order in disk positions appearing for v > v c . 

To investigate the transition in dynamic behavior of grains we conducted 
a set of simulation shearing the static 24x24 configurations, in couette flow. 
Here we show results of simulations with fi = 0.5. (Non-frictional grains and 
other values of friction show similar behavior with the transition occuring at 
v c „, and are thus not presented). The upper wall was moved at a constant 
velocity (v = 10~ 3 £oAo) while / was kept constant, maintaining CVBC. We 
observe different behavior as a function of the solid fraction: In configurations 
with v < u c , momentum is transferred from the wall to interior grains mostly 
via short-lived collisions, observed in spiky fluctuations in stress measured on 
the wall and in coordination number, Fig. [| Stresses are transmitted only 
within local clusters (as in Figure [l|(t 2)). The power spectra of the stress fluc- 
tuations time series approaches white noise demonstrating the uncorrelated 
nature of stress-transfer. For dense packings, with v > u c , grains interact via 
long-lasting contacts. Global motion is characterized by elastic-plastic cycles: 
clusters of grains in contact accumulate recoverable elastic strain, but when 
stresses become too great, grains suddenly rearrange to relieve the stress. 
Continued shearing then begins accumulation of elastic strain on the new 
particle arrangement. This behavior leads to a stick-slip stress time series 
(Figure §). Long clusters of these grains in contact form system-spanning 
'stress chains' that transmit forces from the boundaries into the interior (as 
in Figure |l|(tl)). The stress fluctuation power spectra for v > v c follows f v , 
where rj ~ —2, indicating long-time correlations and in agreement with ex- 
perimental results || conducted in densely packed systems. Most interesting 
is the behavior at the transition point: When v = v c the system oscillates 
between a solid-like 'jammed' state (Figure |l|(tl)) and a gas-like behavior 
(Figure [TJ(t2)). At v — z/ c , values of Z,a,N, fluctuate between values char- 
acteristic of the 'solid-phase' and those characteristic of 'gas-phase' (Figure 
0). The transition density thus bridges the gap between liquid and solid by 
coexistence of the two phases. We note the relation to Kirkwood- Alder tran- 
sition (KAT), a disorder-order transition of repulsive hard-spheres studied 
extensively in the context of collidol suspensions [JR]]. The relation to KAT, 
as well as the transition characteristics, suggest that the granular rigidity 
transition is a first-order transition. 

The relation between stress and strain rate (to be published elsewhere) 
is another property that changes across the phase boundary. When v < v c 
we measure (<r, N) oc v 2 , as expected from theory || and experiments [|TT 



When v > v c measured average stresses are nearly independent of strain rate, 
as expected for elastic-plastic materials and as seen in experiments in densely 
packed granular systems ||12|| . For v — v c stress-strain-rate curves resemble 
those of plastic materials, since high stresses in jammed states dominate 
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time-averaged behavior. 

We finally demonstrate the role of the critical solid fraction when the 
system is allowed to evolve to it's own prefered solid fraction. Soil mech- 
anists have long known of a 'critical density'. If a granular aggregate is 
over-compacted and sheared under CFBC, it will deform while shearing and 
expand to this 'critical density'. If it is initially under-consolidated it will 
compact while shearing until it reaches the same 'critical density'. We per- 
form simulations, similar to these experimental conditions, by applying a 
constant normal stresses to the rigid boundaries (CFBC), while shearing the 
upper wall at velocity v. The systems expand or contract, depending on 
the initial porosity, and finally reach a steady state where variables (u, a, Z) 
fluctuate around a constant value. Steady state solid fractions, < v >, are 
shown in Figure |](a) as function of N and v , for simulations with fi = 0.5. 
For a wide range of applied normal stresses (range widens with decreas- 
ing velocity) systems attain the critical volume fraction for frictional grains, 
i.e. < v >~ u c05 = 0.805. For smooth disks, simulations converge on 
< v >~ u cQ = 0.835, in a similar manner. Though < v > is fairly indepen- 
dent of N in the 'critical regime', the average coordination number < Z > 
increases (Figure f|(bl)(b2)). In those systems having < v >~ v c system- 
spanning stress-chains coexist with unstressed 'gas islands'. With increasing 
N, the number and size of gas islands decreases, and the connectivity of 
chains increases, and thus < Z > increases. Deviations from the critical 
state occur in two cases: 1) At high enough normal stresses, 'solidification' 
occurs: < v > increases above u c , islands virtually vanish and chains become 
fully connected (N = 10~ 2 in Figure |]a). 2) When the ratio of inertial to nor- 
mal forces is high inertial forces cause decompaction, < v > becomes smaller 
than v c (simulations performed at higher velocities and lower normal stresses: 
Figure |]a, v = 10~ 3 and iV < 10~ 5 ), with islands growing to divide stress 
chains, resulting in 'liquification'. We also observed the gas to solid transi- 
tion in the power spectra of stress fluctuations, where power follows f v with r\ 
increasing continuously from r\ ~ for N = 10 -7 to rj ~ —2 for iV > 10 -2 . In 
these CFBC simulations the system, although at the critical state, does not 
flip between gas and solid phases as in CVBC, since 'jamming' episodes may 
be avoided by slight dilation (producing fluctuations around mean porosity 
and height). Instead of temporal coexistence seen in CVBC the critical state 
in CFBC is marked by spatial coexistence of two phases: gas-islands and 
stress-bearing chains. 

Why is the transition between gas and solid also the 'critical-density' to 
which the system is attracted? any finite normal stress acting on grains would 
cause grain compaction and contact, and thus loss of 'gas-like' properties, 
approaching v c from below. Solid fractions of v > v c are characterized by 
finite elastic deformation (> 1%) of grains, which require extremely high 
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normal stresses for stiff natural granular material as rocks. Thus v c is the rigid 
limit, where stresses are accommodated by efficient load bearing structures. 
Experimentally, a coexistence of liquid and solid regimes has been observed 
[[Uj in shearing granular systems, and oscillations between 'jamming' and 
'flowing' states occur spontaneously in a variety of systems, from hoppers 



to natural and experimental land-slides [H, suggesting proximity of these 
systems to the phase boundary. An attracting phase boundary may also 
explain 'fragile materials', a recent term used to describe 'jammed' states 



which may be unjammed by small fluctuations |]T5 



To summarize, we identified the rigidity transition for granular media 
in 2D at v c = 0.80 — 0.84. Two-dimensional results can be mapped to 
corresponding 3D volume fractions using a relation for inter-particle voids, 
v 3D = 4z/^ 2 / 37rl/2 PL predicting that in 3D v c = 0.54 - 0.58 with the 



lower number representing frictional grains. Experiments |17| confirm that 



immense stiffening of rapidly shearing grain aggregates occurs at v c « 0.54. 
The critical solid-fraction is a phase-boundary between gas and solid regimes 
of behavior, where the two phases were observed to coexist in space and/or 
time. The critical state is self-organizing, and corresponds to the 'critical 
density' known in soil mechanics: Sheared aggregates tend to this critical 
volume fraction under a wide range of normal loads. We used normal stresses 
ranging between 10 -2 — f 0~ 7 , which using characteristic young modulus for 
Earth materials (e.g. rock), translates to N = 1 — 10 5 KPa, (corresponding for 
soils to burial depth of 1 — 10 5 m). Based on this we suggest that many natural 
granular deformation processes will occur at a solid fraction which constitutes 
a phase boundary, and neither gas nor elasto-plastic descriptions will fully 
capture these systems behavior. However, such criticality in natural systems 
is easily missed, since spatial and temporal averages of stresses tend to be 
dominated by stress chains, and thus have the mark of solid deformation. It 
is clear that much more theoretical and experimental work needs to be done, 
including investigating much larger systems, and examining the critical state 
more closely. 
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(11) (t2) 



Figure 1: Representative instantaneous grain stresses and configurations dur- 
ing shear (here v = 10~ 3 ), in a CVBC, 24x24 system, with \x = 0.5, at times 
(ti) and (t 2 ). A line is drawn through stressed contacts. The system shown 
here has v = 0.80 ~ u c , and thus exhibits rigid and 'gas-like' behaviors in- 
termittently during shearing: At (ti), it is 'jammed' with system-spanning 
stress chain. At time (t 2 ), the system is 'loose', with only local stress clusters. 
Systems with v > u c always look like the snapshot in (tl), while those with 
v < v c always look like (t2). 
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Figure 2: Results from CVBC simulations of static 12x12 and 24x24 grain 
packings as a function of solid fraction, (a) average coordination number per 
interior grain (b) normal stress exerted on horizontal walls (c) shear modulus 
of the aggregate. Solid curves are theoretical predictions 
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Figure 3: Representative grain a) coordination numbers b) shear stresses, as 
a function of strain, for a CVBC 24x24 system, with /i = 0.5, sheared at 
v = 10 -3 . Different lines represent time-series for runs having v < i> c , v = 
v c = 0.80, and v > u c . 
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Figure 4: (a) Time averaged solid fraction of 24x24 systems, with ji = 0.5, 
sheared at two different velocities, under CFBC. Error bars depict standard 
deviation from steady-state, (b) Grain configurations and stresses for two of 
the v = 10~ 4 runs in (a): (bl) applied normal stress: N = 10~ 6 , measured 
time-averaged solid fraction < v >= 0.806 ±0.002 and coordination number: 
< Z >= 1.92 ± 0.13, (b2) N = 10' 3 , < v >= 0.809 ± 0.003, < Z >= 
2.86 ±0.06. 
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